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Anthony  V.  Fiacco 

1.  Introduction 

Initial  numerical  results  resulting  from  the  implementation  of  a 
penalty  function  technique  for  obtaining  sensitivity  information  in  para- 
metric nonlinear  programming  were  given  by  Armacost  and  Fiacco  (1974). 

The  work  is  based  on  the  theory  developed  by  Fiacco  and  McCormick  (1968) 
and  extended  by  Fiacco  (1973).  This  paper  reports  on  refinements  and  ex- 
tensions of  the  computational  procedures  implemented  by  Armacost  and 
Mylander  (1973)  and  Armacost  (1976),  using  the  SUMT-Version  4 computer 
code  with  the  logarithmic-quadratic  loss  penalty  function  to  estimate 
the  partial  derivatives  of  the  solution  point  and  the  objective  function 
optimal  value,  the  derivatives  here  being  taken  with  respect  to  the 
specified  problem  parameters. 

Fiacco  (1973)  developed  the  necessary  general  formulas  for  the 
partial  derivatives  of  the  "optimal  value  function,"  the  components  of  a 
local  solution  point  and  its  associated  optimal  Lagrange  multipliers,  for 
a large  class  of  parametric  nonlinear  programming  problems  composed  of 
twice  differentiable  functions.  He  also  obtained  approximation  formulas 
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in  terms  of  the  well  known  logarithmic-quadratic  penalty  function.  Re- 
cently, Armacost  and  Fiacco  (1975)  particularized  and  simplified  these 
formulas  for  various  problem  structures  and  developed  formulas  for  the  first 
and  second  derivatives  of  the  optimal  value  function  of  the  given  problem. 
Additionally,  Armacost  and  Fiacco  (1976)  have  applied  the  general  theory 
to  easily  prove  the  well  known  result  that,  when  the  parameters  are  the 
right-hand  side  components  of  the  constraints,  the  optimal  Lagrange  multi- 
pliers give  the  gradient  of  the  optimal  value  function  (with  respect  to 
the  parameters).  Further,  it  was  shown  that  the  first  derivatives  of  the 
Lagrange  multipliers  give  the  components  of  the  Hessian  of  the  optimal 
value  function,  and  explicit  formulas  were  developed  for  the  Hessian  in  terms 
of  the  problem  functions. 

In  their  first  report  on  computational  experience,  Armacost  and 
Fiacco  (1974)  concentrated  primarily  on  presenting  computational  experience 
associated  with  the  calculation  of  the  first  derivatives  of  a local  solu- 
tion point.  The  practical  implementability  of  the  approach  was  demonstrated. 

In  a subsequent  paper,  Armacost  (1976)  reported  on  additional  compu- 
tational experience,  focusing  on  the  calculation  of  the  derivatives  of  the 
optimal  value  function  and  the  Lagrange  multipliers,  also  implementing  a 
potentially  valuable  refinement  that  allows  for  computerized  screening  for 
"key"  parameters.  This  paper  may  be  regarded  as  a continuation  and  ampli- 
fication of  the  Armacost  paper. 

For  problems  involving  a large  number  of  parameters,  a very  large 
number  of  partial  derivatives  may  be  calculated  if  one  proceeds  indiscrim- 
inately. This  is  not  only  time-consuming,  but  may  also  be  quite  burdensome 
to  a user  who  must  evaluate  the  overall  significance  of  the  results.  One 
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measure  of  the  latter  is  the  effect  of  a perturbation  on  the  solution  value. 
It  is  quite  possible  and  often  observed  in  practice  that  the  optimal  objec- 
tive function  value  is  much  more  sensitive  to  a few  of  the  many  parameters 
present.  With  this  in  mind,  the  method  developed  by  Armacost  and  Fiacco 
(1975)  to  estimate  the  first  order  sensitivity  of  the  optimal  value  func- 
tion was  incorporated  in  the  computer  program  to  provide  an  option  for 
preliminary  screening  of  the  parameters  to  eliminate  further  calculations 
involving  perturbations  of  parameters  having  "little"  effect  on  the  optimal 
value  function.  (A  user  can  easily  introduce  his  own  criteria  of  signifi- 
cance in  this  determination.)  Using  the  formulas  developed  by  Fiacco  (1973), 
a second  option  is  included  which  permits  the  calculation  of  the  sensitivity 
estimates  for  the  Lagrange  multipliers.  The  computer  code  and  options 
used  to  accomplish  these  and  other  calculations  are  discussed  in  Section  3. 

In  Section  4,  a sensitivity  analysis  is  conducted  for  a multi-item 
inventory  model  developed  by  Schrady  and  Choe  (1971)  for  the  U.S.  Navy. 

The  example  analyzed  is  the  same  small  one  treated  by  Schrady  and  Choe, 
though  readily  extended  to  a large-scale  model.  The  results  illustrate 
the  potential  value  of  a detailed  automated  post-optimality  sensitivity 
analysis  in  practical  situations,  and  hopefully  dramatize  the  numerous 
rich  interpretations  and  insights  that  can  be  derived  from  this  information, 
as  well  as  indicating  the  caution  that  must  be  taken  in  making  valid  in- 
ferences. 

The  recently  obtained  basic  theoretical  results  validating  the  compu- 
tational algorithm  are  summarized  rather  completely  in  the  next  section 
so  that  the  paper  might  be  self  contained. 

- 3 - 
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2.  Supporting  Theory 

The  parametric  mathematical  programming  problems  considered  here 
are  of  the  form 


minimize  f(x,e) 
x e En 


subject  to  g^ (x , £ ) ^ 0 , 
hj (x,e)  = 0 , 


where  x is  the  usual  vector  of  variables  and  £ is  a k-component  vector 
of  numbers  called  "parameters."  It  is  desired  ultimately  to  develop  a 
complete  characterization  of  a solution  x(e)  of  Problem  P(e)  as  a function 
of  £ . In  our  current  work,  we  have  concentrated  on  certain  recently 
computationally  tractable  measures  of  change  in  a solution  as  £ is  per- 
turbed from  a specified  value.  (Without  loss  of  generality,  we  assume 
that  the  specified  value  is  £=0  .) 

When  certain  assumptions  are  satisfied,  Fiacco  (1973)  and  Armacost 
and  Fiacco  (1975)  have  characterized  the  "first  order  sensitivity"  of  a 
"Kuhn-Tucker  Triple"  and  the  first  and  second  order  sensitivity  of  the 
optimal  value  function  of  Problem  P(e).  (These  quantities  are  defined 
as  the  theory  is  presented.)  Additionally,  they  have  developed  formulas 
for  efficiently  estimating  this  sensitivity  when  the  logarithmic-quadratic 
loss  penalty  function  algorithm  is  used  to  solve  Problem  P(e).  The  main 
theoretical  results  are  summarized  here. 

The  Lagrangian  for  Problem  P(£)  is  defined  as 


m p 

L(x,u,w,e)  = f(x,e)  - l u g (x,e)  + \ w.h  (x,e)  , 

i=l  1 1=1  J J 
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where  ik  , and  , j=l,...,p  are  "Lagrange  multipliers 

associated  with  the  inequality  and  equality  constraints,  respectively. 

Any  vector  (x,u,w)  satisfying  the  usual  (first  order)  Kuhn-Tucker  condi- 
tions (Fiacco  and  McCormick,  1968)  of  Problem  P(e)  is  called  a Kuhn-Tucker 
triple. 

The  following  four  assumptions  are  sufficient  to  establish  the 
results  and  are  assumed  to  hold  throughout  the  paper: 

A1  — The  functions  defining  Problem  P(e)  are  twice  continuously 
differentiable  in  (x,e)  in  a neighborhood  of  (x*,0)  . 

A2  — The  second  order  sufficient  conditions  for  a local  minimum 
of  Problem  P(0)  hold  at  x*  with  associated  Lagrange  mul- 
tipliers u*  and  w*  . 

A3  — The  gradients  Vxg^(x*,0)  (i.e.,  (9g^(x*,0) /3x^ , . . . , 

T 

9g^ (x*, 0)/3xn>  , the  superscript  T denoting  transposition) 

for  all  i such  that  g. (x*,0)  = 0 , and  V h.(x*,0)  , 

i x 3 

j=l,...,p  are  linearly  independent. 

A4  — Strict  complementary  slackness  holds  at  (x*,0)  (i.e., 

u*  > 0 for  all  i such  that  g^(x*,0)  = 0 ). 

Theorem  1:  (Local  characterization  of  a Kuhn-Tucker  triple  (Fiacco, 

1973)  of  Problem  P(e).)  If  assumptions  Al,  A2 , A3  and  A4  hold  for  Problem 
P(e)  at  (x*,0)  , then 

(a)  x*  is  a local  isolated  minimizing  point  of  Problem  P(0) 
and  the  associated  Lagrange  multipliers  u*  and  w*  are 
unique ; 

(b)  for  e in  a neighborhood  of  0 , there  exists  a unique,  once 
continuously  differentiable  vector  function  y(£)  - 
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T 

(x(e) ,u(e) ,w(e))  satisfying  the  second  order  sufficient 

conditions  for  a local  minimum  of  Problem  P(e)  such  that 
T 

y(0)  = (x*,u*,w*)  = y*  and  hence,  x(e)  is  a locally 

unique,  local  minimum  of  Problem  P(e)  with  associated  unique 
Lagrange  multipliers  u(e)  and  w(e)  ; and 
(c)  for  £ near  0 , the  set  of  binding  inequalities  is  unchanged, 
strict  complementary  slackness  holds  for  u^(e)  for  i such 
that  g^(x(e),£)  = 0 , and  the  binding  constraint  gradients 
are  linearly  independent  at  x(e)  . 

This  result  provides  a characterization  of  a local  solution  of 
Problem  P(e)  and  its  associated  optimal  Lagrange  multipliers  near  £=0  . 

It  generalizes  a theorem  first  presented  by  Fiacco  and  McCormick  (1968, 
Theorem  6)  and  is  closely  related  to  a generalization  of  the  same  theorem 
provided  independently  by  Robinson  (1974).  It  shows  that  the  Kuhn-Tucker 
triple  y(e)  is  unique  and  well  behaved,  under  the  given  conditions. 

Since  y(e)  is  once  differentiable,  the  partial  derivatives  of  the  compo- 
nents of  y(e)  are  well  defined.  This  fact  and  Assumption  A1  also  mean 
that  the  functions  defining  Problem  P(e)  are  once  continuously  differentiable 
functions  of  £ along  the  "solution  trajectory"  x(e)  near  £=0  , and 
the  Lagrangian  is  a once  continuously  differentiable  function  of  £ along 
the  "Kuhn-Tucker  point  trajectory." 

We  are  thus  motivated  to  determine  a means  to  calculate  the  various 
partial  derivatives,  since  this  yields  a first  order  estimate  of  the  locally 
optimal  Kuhn-Tucker  triple  and  the  problem  functions  near  £=0  . 

Denote  by  V x(£)  = (9x.(e)/9e.)  , i=l,...,n  , j=l,...,k  , the 
® i J 

n x k matrix  of  partial  derivatives  of  x(e)  with  respect  to  £ , and 
define  V^u(e)  and  V^w(e)  in  a similar  fashion.  We  then  define 

- 6 - 
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Vey(e)  = (V^x(e) ,V^u(e) ,V^w(e))T  , an  (n+m+p)  x k matrix. 

When  y(e)  is  available,  V^y(e)  can  be  calculated  by  noting  that 
Conclusion  (b)  of  the  theorem  implies  the  satisfaction  of  the  Kuhn-Tucker 
conditions  for  P(e)  at  y(e)  near  e=0  , i.e. , 

VxL[x(e) ,u(e) ,w(e) ,e]  = 0 , 

ui(e)g1[x(e) ,£]  = 0 , i=l,...,m  , (1) 

h^[x(c),£]  — 0 , j— l,...,p  • 

Since  the  Jacobian  M(e)  of  this  system  with  respect  to  (x,u  ,w)  (i.e., 

the  matrix  obtained  by  differentiating  the  left  side  of  (1)  with  respect 
to  the  components  of  (x,u,w)  ) is  nonsingular  under  the  given  assumptions, 
the  total  derivative  of  the  system  with  respect  to  e is  well  defined  and 
must  equal  zero.  This  yields 

M(e)Vcy(e)  = N(e)  , 

where  N(e)  is  the  negative  of  the  Jacobian  of  the  Kuhn-Tucker  system 
with  respect  to  z , and  hence 

Vcy(e)  = M(e)  LN(e)  . 

The  class  of  algorithms  based  on  twice  continuously  differentiable 
penalty  functions  can  be  used  without  additional  assumptions  and  without 
requiring  y (e)  to  provide  an  estimate  of  V^yCe)  . Furthermore,  most  of 
the  information  required  to  make  the  estimate  is^  already  available  in  the 
typical  implementations  of  these  algorithms.  Here,  we  use  the  logarithmic- 
quadratic  penalty  function  for  Problem  P(e)  (Fiacco  and  McCormick,  1968) 
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iu  r „ 

(x,e, r)  = f(x,e)  - r £ £n  g (x,c)  + (l/2r)  £ h.(x,e)  . (2) 

i=l  j=l  J 


Under  the  given  assumptions,  the  following  facts  are  known  for 
Problem  P(0)  from  penalty  function  theory  (Fiacco  and  McCormick,  1968, 
Theorems  10  and  17): 

(1)  For  r > 0 and  small,  there  exists  a unique  once  continu- 
ously differentiable  vector  function  x(0,r)  such  that 
x(0,r)  is  a locally  unique  minimizing  point  of  W(x,0,r) 
in  R^(0)  = {x:  g^(x,0)  > 0 , i=l,...,m  , and  h.(x,0)  = 0 , 
j=l,...,p}  and  such  that  x(0,r)  -*■  x(0,0)  = x*  ; 


(2)  lim  r l in  g [x(0,r)]  = 0 ; 
r-*0  i=l 


(3)  lim  (l/2r)  J h.[x(0,r),0]  = 0 ; and 


(4)  lim  W[x(0,r) ,0,r ] = f(x*,0) 
r-K) 


The  following  theorem  extends  these  results  for  Problem  P(£),  where 

£ is  allowed  to  vary  in  a neighborhood  of  0 , and  provides  a basis  for 

approximating  the  sensitivity  information  associated  with  Problem  P(e). 

2 

The  notation  V W denotes  the  matrix  of  second  partial  derivatives  of  W 
x r 

with  respect  to  x . 

Theorem  2 : (Relationship  of  solutions  of  Problem  P(e)  and  minima 

of  W(x,e,r)  , (Fiacco  1973).)  If  Assumptions  A1  - A4  hold,  then  in  a 

neighborhood  about  (£,r)  = (0,0)  there  exists  a unique  once  continuously 

T 

differentiable  vector  function  y(£,r)  = [x(£ ,r) ,u(£ ,r)w(e ,r) ] satisfying 
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VxL(x,u,w,e)  = 0 , 

uigi(x,e)  = r , i=l , . . . ,m  , 

hj (x,e)  = w^r  , j=l,. . . ,p  , 

with  y(0,0)  = (x*,u*,w*)  and  such  that,  for  any  (e,r)  near  (0,0) 

and  r > 0 , x(e,r)  is  a locally  unique  unconstrained  local  minimizing 

2 

point  of  W(x,e,r)  , gi[x(e,r),e]  > 0 , i=l,...,m  , and  VxW[x(e,r) ,e,r] 
is  positive  definite. 

Corollary  2.1:  (Convergence  of  estimates  using  W(x,e,r)  , (Fiacco, 

1973).)  If  Assumptions  Al,  A2,  A3  and  A4  hold  for  Problem  P(e) , then  for 
any  £ near  0 , 

(a)  lim+  y(e,r)  = y(e,0)  = y(e)  , the  Kuhn-Tucker  triple  charac- 
r-K) 

terized  in  Theorem  1;  and 

(b)  lim  V y(e,r)  = V y(e,0)  = V y(e)  . 

r-*0  K e 

This  result  motivates  use  of  V^y(e,r)  to  estimate  VFy(e)  , when 
e is  near  0 and  r is  near  0 , once  y(e,r)  is  available.  Theorem  2 
provides  the  basis  for  an  efficient  calculation  of  V£y(e,r)  . Since,  at 
a local  solution  point  x(e,r)  of  W(x,e,r)  , it  follows  that 

VxW[x(e,r) ,e,r]  = 0 , (3) 

we  can  differentiate  (3)  with  respect  to  e to  obtain 

V^W[x(e,r) ,e,r]V  x(e,r)  + V (V  W[x(e,r) ,e,r])  = 0 . (4) 

X c.  EX 

2 

By  Theorem  2,  V^W  is  positive  definite  for  (e,r)  near  (0,0) 

2 2-1 
and  r > 0 , so  V W has  an  inverse  and  V x(e,r)  = -V  W[x(e,r) ,e,r] 

X G X 


luimmwwip  i .iiw'im 


• . ■ 


T-340 


Also,  since 

ut(e,r)  = r/gi(x(e,r) ,E)  , i=l,...,m  , (5) 

and 

w^(e,r)  = (x(r.,r)  ,e)/r  , j = l,...,p  , (6) 

for  (e, r)  near  (0,0)  and  r > 0 , these  equations  can  be  differentiated 
with  respect  to  e to  obtain 

V£.ui(e,r)  = -(r/g^)  [Vxg1(x(e,r)  ,e)  • V^xCe.r)  + dgi(x(e,r)  ,e)/9e]  , 

(7) 

V w.(e,r)  = (1/ r) [V  h.(x(e,r),E)  • V x(E,r)  + 9h. (x(e,r) ,e)/9e]  . (8) 

£ J x 1 £ 3 

Solving  (4)  and  calculating  (7)  and  (8)  then  yields  the  components  of 
V£y(e,r)  , which  can  be  used  to  estimate  V^y (e)  for  (e,r)  near  (0,0)  . 

The  next  results  extend  this  theory  to  an  analysis  of  the  optimal 
value  function  of  Problem  P(e)  along  the  Kuhn-Tucker  point  trajectory 
[x(e) ,u(e) ,w(e)]T  . 

The  optimal  value  function  is  defined  as: 

f*(£)  = f[x(E),e]  , (9) 

and  the  "optimal  value  Lagrangian"  is  defined  as: 

L*(e)  = L[x(e),u(e),w(e),e]  . (30) 

Theorem  3:  (First  and  second  order  changes  in  the  optimal  value 

function,  Armacost  and  Fiacco  (1975).)  If  assumptions  A1  - A4  hold  for 
Problem  P(e),  then  for  £ near  0 , f*(E)  is  a twice  continuously 

differentiable  function  of  £ , and 
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(a)  f*(e)  = L*(e)  ; 


(b)  V f *(e)  = V L(x,u,w,e) 

(x,u,w)=(x(e) ,u(e) ,w(e)) 
ra 

= V f(x,£)  - l u V g (x,e) 
e i=l  1 e i 

P 

+ £ w. V h . (x,e) 

j=l  3 e 3 (x,u,w)=(x(e),u(e),w(e))  ; 

(c)  V2f*(e)  = Ve(VeL(x(e),u(e),w(e),e)T)  . 

The  logarithmic-quadratic  loss  penalty  function  (2)  can  also  be 
used  to  provide  estimates  of  the  first  and  second  order  sensitivity  of 
the  optimal  value  function.  Let  the  optimal  value  penalty  function  be 
defined  as  W*(e,r)  = W(x(£,r) ,£,r)  . 

Theorem  4 : (First  and  second  order  sensitivity  of  W*(e,r)  and 

estimates  for  f*(£)  , Armacost  and  Fiacco  (1975).)  If  Assumptions 
A1  - A4  hold  for  Problem  P(e),  then  for  (e,r)  near  (0,0)  and  r > 0 , 
W*(e,r)  is  a twice  continuously  differentiable  function  of  £ and 

(a)  lim  W*(£ ,r)  = L*(e)  = f*(£)  ; 
r->0 

(b)  V W*(£,r)  = V L(x,u,w,e) 

(x,u,w)=(x(e,r) ,u(E,r) ,w(e,r))  ; 

(c)  lim,  V W*(e , r)  = V L(x(e) ,u(e) ,w(e) ,e)  = V f*(£)  ; (11) 

r-0+  £ £ £ 

(d)  W*(£,r)  = Ve(V£L(x(£,r) ,u(£,r) ,w(e,r) ,e)T)  ; 

(e)  lim  V2W*(£,r)  = V2f*(£)  . 
r+0  E 

This  result  provides  a justification  for  estimating  f*(£)  , 

V£f *(£)  and  V2f*(£)  by  W*(e,r)  , VfW*(£,r)  and  V2W*(£,r)  , respectively 


when  r is  positive  and  small  enough. 
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Since  Corollary  2.1  and  continuity  imply  that  lim+  f(x(e,r),e)  = 

r-K) 

f*(e)  , another  estimate  of  the  optimal  value  function  (9)  is  provided  by 
f'(e,r)  = f (x(e,r) ,e)  when  r > 0 and  small.  Direct  application  of  the 
chain  rule  for  differentiation  then  yields,  for  x = x(c, r)  , 


(12) 


Under  the  given  assumptions,  continuity  also  assures  that  V^f  (e,r)  -*■ 
V£f*(c)  as  r -*•  0+  . Thus,  both  V£f^(e,r)  and  Vf,W*(£,r)  are  estimates 

of  V f*(e)  for  r sufficiently  small. 

e 

It  should  be  noted  that  these  estimates  are  functionally  related 

since 


m 


V W*(e,r)  = V f'(e,r)  - l u± (Vxg±V  x(e , r)  + V^g.) 

i=l 


+ l w.(V  h.V  x(e,r)  + V h.) 

j“1  J * j e e j 


x=x(e,r) 


« 


i 


j 


l 1 


t 


From  this  expression,  it  is  clear  that  V^f  (e,r)  is  the  better  estimate 
of  V f*(e)  , the  remaining  terms  in  V^W*(e,r)  simply  constituting  "noise" 
that  is  eliminated  as  r -*■  0+  . However,  by  using  the  expression  for 
V^W*(£,r)  given  by  (11),  V^W*(e,r)  can  be  evaluated  without  necessitating 

the  calculation  of  Vex(e,r)  , which  is  required  to  compute  (12).  Thus, 
the  cruder  but  computationally  much  cheaper  estimate  of  V£f*(£)  given  by 
Equation  (11)  has  now  been  introduced  as  an  option  in  the  computer  program 
as  a preliminary  screening  device  to  identify  crucial  parameters.  Restric- 
tion of  subsequent  calculations  to  these  parameters,  and  other  calculations 
such  as  the  sharper  estimate  of  f*(e)  given  by  (12)  are  provided  as 
additional  options. 
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In  summary,  the  basis  for  the  estimation  procedure  utilized  here 
for  a specific  problem,  say  Problem  P(0),  is  the  minimization  of  the  penalty 
function  W(x,e,r)  given  by  (2).  This  yields  a point  x(e,r)  which 
may  be  viewed  as  an  estimate  of  a (local)  solution  x*  of  Problem  P(0). 

The  associated  optimal  Lagrange  multipliers  u*  and  w*  are  then  estimated 
by  using  the  relationships  given  in  (5)  and  (6),  respectively.  The  first 
partial  derivatives  of  these  quantities  with  respect  to  e are  then 
obtained  by  first  solving  (4)  and  then  applying  (7)  and  (8).  The  estimate 
f(x(e,r),e)  of  f*(e)  is  immediately  available  when  x(e,r)  has  been 
determined,  and  the  two  estimates  of  V^f*(e)  given  by  (11)  and  (12) 
were  already  discussed.  Various  options  and  computer  codes  implementing 
the  procedure  are  discussed  in  the  next  section. 

3.  User  Options,  Computer  Codes  and  an  Example 

The  basic  SUMT-Version  4 computer  program  and  instructions  for 
its  use  are  described  in  Mylander,  Holmes  and  McCormick  (1971).  The 
basic  sensitivity  analysis  subroutines,  user  instructions,  and  instruc- 
tions for  integrating  the  sensitivity  package  with  the  SUMT-Version  4 
code  are  described  in  Armacost  and  Mylander  (1973). 

Briefly,  the  conduct  of  a sensitivity  analysis  is  controlled  by  the 
variable  NEX0P3  which  is  given  a value  on  the  "Second  Option  Card"  in  the 
SUMT  input  data  deck.  There  are  four  choices:  no  sensitivity  analysis, 

a sensitivity  analysis  at  the  final  subproblem,  a sensitivity  analysis  at 
each  subproblem  along  the  penalty  function  minimizing  trajectory,  or  a 
sensitivity  analysis  at  the  final  subproblem  for  a range  of  differencing 
increments.  In  conjunction  with  this  option,  two  additional  options  are 
added  here  and  come  into  play  whenever  a sensitivity  analysis  is  conducted. 
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The  first  option  (technically  Option  4)  is  controlled  by  the 
variable  NEX0P4  and  determines  whether  the  partial  derivatives  of  the 
Lagrange  multipliers  will  be  calculated.  When  the  calculation  is  done, 
the  formulas  developed  by  Fiacco  (1973)  are  used. 

The  second  option  added  here  (Option  5)  permits  a screening  of  the 
parameters  to  reduce  the  number  of  partial  derivatives  which  are  estimated 
by  limiting  further  analysis  to  parameter  changes  which  affect  the  opti- 
mal value  of  the  objective  function  by  an  amount  exceeding  0.1  percent  of 
its  current  value.  This  option  is  controlled  by  the  variable  NEX0P5. 

The  estimate  of  sensitivity  of  the  optimal  value  function  with  respect  to 
a particular  parameter  under  this  option  is  calculated  using  the  Armacost 
and  Fiacco  (1975)  result  which  involves  the  partial  derivative  of  the 
Lagrangian  taken  with  respect  to  the  parameter  under  consideration. 

Subroutines  LMULT  and  PRESEN  and  related  coding  in  Subroutine  SENS 
implement  Option  4 and  Option  5,  respectively.  Subroutines  SENS,  LMULT 
and  PRESEN  are  listed  in  Appendix  A.  Specific  instructions  for  using 
these  two  options  in  conjunction  with  the  "Second  Option  Card"  are  given 
below  in  Table  1.  This  information  should  be  added  to  Table  5 in  My lander. 
Holmes  and  McCormick  (1971). 

As  an  illustration  of  the  kind  of  information  that  can  be  generated, 
consider  the  following  simple  (convex)  parametric  nonlinear  programming 
problem. 


minimize 

f(x,e)  = Xj^  + e2x2 

subject  to 

, v 2 2 

81(x,e)  * Ej  - x • 

for  e1  > 0 and  e2  not  restricted. 


A(e) 
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TABLE  1 

THE  SECOND  OPTION  CARD 


Option 

Column 

Value 

Meaning 

4 

28 

=0 

Do  not  estimate  the  partial  derivatives  of  the 
estimates  of  the  Lagrange  multipliers. 

=1 

Estimate  the  partial  derivatives  of  the  esti- 
mates of  the  Lagrange  multipliers  whenever  a 
sensitivity  analysis  of  the  solution  point  is 
conducted. 

5 

35 

=0 

Estimate  the  partial  derivatives  of  the  opti- 
mal value  function  and  eliminate  those  param- 
eters which  do  not  affect  the  optimal  value 
function  from  subsequent  sensitivity  calcula- 
tions. 

=1 

Estimate  the  partial  derivatives  of  the  opti- 
mal value  function  with  respect  to  all  param- 
eters, but  continue  all  subsequent  sensitivity 
calculations  with  respect  to  all  parameters. 

=2 

Do  not  estimate  the  partial  derivatives  of  the 
optimal  value  function  first.  Conduct  the 
sensitivity  analysis  with  respect  to  all  pa- 
rameters . 

To  dramatize  the  complexity  of  the  relationships  that  can  arise 
even  in  such  trivial  problems,  we  first  give  the  solution  in  closed  form. 
Application  of  the  Kuhn-Tucker  conditions  yields  the  following: 


x(e) 


x-^e) 

x2(e) 


1 + et 


he2 


V 


1 + e; 
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u(e)  = ■yl  + e|  / 2 e1  , 


and 


g1(x(e) ,e)  = 0 . 


(It  may  be  readily  verified  that  Assumptions  A1  - A4  are  satisfied.) 
From  the  expression  for  x(e)  , we  find  that 


f*(e)  = -EjY1  + e2  . 


Taking  partial  derivatives  with  respect  to  e and  , we  obtain 


V x(e)  = 

e 


+ e; 


ele2 


(1  + e ) 


-e. 


(1  + e?) 


V f*(e)  = (3f*(E)/3e1,  3f*(e)/3e2) 


-yi  + e2  , -e1e2/ %/ 1 + e2  I * 


and 


V£u*(e)  = 


V1  + 4 nH  • e2/(2elV1  + e2  )) 


Suppose  we  are  particularly  interested  in  the  values  e^=2  » e2=1 
Evaluation  of  the  above  expressions  yields 


f*  = -iSi  , V^f*  = (-/I,  -Si)  , 


~-Si 

-Sin 

Sin 

X*  = 

-Si_ 

, V£x*  = 

-Sin 

-Sin 

jk.; 

• 

- 16  - 

T-340 


u*  = Jl!  4 


V£u * = (-^2/8,  ^2/8)  . 


The  numerical  results  obtained  by  the  computer  program  are  included 
in  Table  2 for  the  optimal  value  function  and  Lagrange  multiplier  sensi- 
tivity. The  values  of  the  first  order  optimal  value  function  sensitivity 
computed  both  by  the  chain  rule  (Equation  (12))  and  by  taking  partial 
derivatives  of  the  Lagrangian  with  respect  to  the  parameters  (Equation  (11)) 
are  presented  for  comparison.  As  indicated  previously  and  expected,  the 
chain  rule  estimates  are  sharper  The  subproblems  correspond  to  minimiza- 
tions of  the  penalty  function  W(x,e,r)  for  several  values  of  r > 0 , 
the  value  of  r decreasing  as  the  subproblem  index  increases.  The  values 
of  f are  given  by  f(x(e,r),e)  , while  the  values  of  u and  V^u  were 
obtained  by  applying  formulas  (5)  and  (7).  The  components  of  x(e)  and 

V x(e)  were  calculated  also  but  are  not  portrayed. 

E 


TABLE  2 

SENSITIVITY  ESTIMATES  FOR  PROBLEM  A(e) , WHERE  e=(2,l) 


Subproblem 

f 

Lagrangian 
9f/9s1  9f/9e2 

Chain 

9f/9e1 

Rule 

3f/9e2 

u 

9u/9e^ 

9u/ 9e2 

1 

-1.9999 

-1.9999 

- .9999 

-1. 3333 

-1. 3333 

.4999 

-.3333 

.1666 

2 

-2.3393 

-1.5440 

-1.2947 

-1.4087 

-1.4088 

.3860 

-.2100 

.1761 

3 

-2.7765 

-1.4439 

-1.3833 

-1.4139 

-1.4139 

.3609 

-.1844 

.1767 

4 

-2.8128 

-1.4243 

-1.4064 

-1.4142 

-1.4142 

. 3560 

-.1790 

.1768 

5 

-2.8245 

-1.4127 

-1.4123 

-1.4142 

-1.4142 

.3531 

-.1768 

.1768 

6 

-2.8274 

-1.4137 

-1.4137 

-1.4142 

-1.4142 

. 3532 

-.1767 

.1768 

7 

-2.8282 

-1.3899 

-1.4141 

-1.4142 

-1.4142 

.3475 

-.1737 

.1768 

Analytical 

-2.8282 

-1.4142 

-1.4142 

-1.4142 

-1.4142 

.3537 

-.1769 

.1768 
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The  results  are  typical  of  those  that  have  been  obtained  to  date, 
convergence  of  the  sensitivity  estimates  being  reliable  and  stable  with 
less  stability  noted  in  the  estimates  of  the  optimal  Lagrange  multipliers. 
It  appears  that  the  sensitivity  estimates  will  converge  satisfactorily 
whenever  the  SUMT  algorithm  succeeds  in  determining  a good  estimate  of  a 
Kuhn-Tucker  point. 

4.  A Large-scale  Multi-item  Inventory  Model 

Traditionally,  inventory  models  have  been  formulated  to  minimize 
some  function  of  the  ordering,  holding  and  shortage  (or  backorder)  costs 
subject  to  various  constraints.  Schrady  and  Choe  (1971)  have  formulated 
an  inventory  model  which  appears  to  have  much  greater  relevance  for  an 
inventory  system  in  a noncommercial  environment,  such  as  institutional 
or  military.  The  costs  used  in  the  traditional  models  may  be  quite  arti- 
ficial and  the  real  objective  of  the  system  is  often  maximization  of  a 
measure  of  readiness  or  service,  here  assumed  to  be  equivalent  to  minimiza- 
tion of  stockouts.  In  addition,  the  stock  points  of  such  supply  systems 
are  inevitably  constrained  by  investment  and  reorder  work  load  limitations. 

Schrady  and  Choe's  multi-item  inventory  system  assumes  these  con- 
straints along  with  the  specific  objective  of  minimizing  the  total  time- 
weighted  shortages.  The  decision  variables  are  taken  to  be  the  "reorder 
quantities"  and  the  "reorder  points,"  respectively,  how  much  to  order  and 
when  to  order  each  item  in  the  inventory.  A three-item  example  problem 
was  solved  by  Schrady  and  Choe  (1971)  using  the  SUMT  computer  code 
(Mylander,  et^  a^.  , 1971).  Subsequently,  McCormick  (1972)  showed  how  the 
special  structure  of  this  inventory  model  can  be  used  to  facilitate  the 
use  of  the  SUMT  code  to  solve  very  large  inventory  problems.  He  also 
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extended  the  model  to  include  constraints  on  storage  volume  and  the  prob- 
ability of  depletion  of  critical  items. 

The  model  and  example  presented  here  are  the  original  ones  due  to 

V 

Schrady  and  Choe.  The  computer  program  described  in  the  preceding  section 
was  used  to  solve  the  example  and  calculate  the  partial  derivatives  of 
various  quantities  of  interest,  with  respect  to  each  parameter  involved 
in  defining  the  model.  (The  analysis  can  be  applied  to  the  extended  model 
without  difficulty. ) 

Detailed  development  of  the  model  is  beyond  the  scope  of  this  paper. 
The  interested  reader  is  referred  to  the  Schrady-Choe  and  McCormick  papers. 
Here,  we  give  a summary  treatment  of  the  various  conditions  and  relation- 
ships upon  which  the  model  is  based.  We  then  tabulate  the  results  obtained 
in  solving  the  resulting  nonlinear  programming  problem  and  applying 
the  sensitivity  analysis  methodology.  A number  of  observations  and  in- 
terpretations are  offered  to  illustrate  the  many  uses  to  which  the  sensi- 
tivity information  might  be  applied. 

It  is  assumed  that  the  amount  of  each  item  in  inventory  is  always 
known,  that  all  demand  which  occurs  when  the  on-hand  stock  is  zero  is 
backordered,  and  that  the  demand  which  occurs  during  the  time  between  the 
placement  of  an  order  and  its  receipt  by  the  stock  point  (i.e.,  the  "lead 

2 

time  demand")  is  normally  distributed  with  known  mean  and  variance  . 

For  the  ith  item,  let 

c^  = item  unit  cost  (in  dollars) , 

A^  = mean  demand  per  unit  time  (in  units), 
r^  = reorder  point, 

= reorder  quantity, 
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The  expected  on-hand  inventory  of  item  i is  given  by  r + Q^/2  - + 

and  the  expected  number  of  orders  placed  per  unit  time  for  item 

i is  A./Q.  . 
r i 

Using  the  above  expressions  and  assumptions,  Schrady  and  Choe  (1971) 
indicate  that  meaningful  approximations  of  the  given  quantitites  are 
obtained  even  when  the  second  term  is  dropped  from  the  expression  for  the 
expected  shortages,  and  when  the  last  term  is  dropped  from  the  expression 
for  expected  on-hand  inventory.  The  given  assumptions  and  simplifications 
then  lead  readily  to  the  following  nonlinear  programming  problem  (which 
Schrady  and  Choe  1971)  proved  convex), 


Wmm 


(SC) 


N 

subject  to  g^Q.r)  = - £ c.^  + Qj/2  " l^)  > 0 , 

i=l 

N 

g2(Q,r)  = K2  - l X./Q.  > 0 , 
i=l 

T 

with  unrestricted  in  sign,  0 , i=l,...,N  , Q=  (Q^,...,QN)  , 

T 

r = (r^,...,rN)  , and  g^  and  g2  representing  the  investment  and 

work  load  constraints,  respectively. 

The  problem  data  for  the  Schrady-Choe  three-item  example  and  the 
initial  starting  point  for  the  SUNT  program  are  shown  in  Table  3.  As 
indicated  in  the  table,  the  lead  time  demands  and  standard  deviations, 
the  item  unit  costs  and  mean  demands,  and  the  investment  and  work  load 
limits  are  all  treated  as  parameters  in  conducting  the  sensitivity  analysis. 

Table  4 gives  the  computer  solution  and  Table  5 the  final  estimate 
of  the  first  partial  derivatives  of  the  optimal  value  function  Z*  with 
respect  to  the  problem  parameters.  Relative  to  the  criterion  used  in  the 
computer  program,  the  Table  5 results  indicate  that  the  optimal  value 
function  is  sensitive  to  parameters  , O and  c^  . 

Many  inferences  are  possible.  For  example,  the  fact  that  the  solution  is 
particularly  sensitive  to  the  values  of  the  standard  deviations  of  the 
lead  time  demand  of  items  2 and  3 might  indicate  that,  since  these  parameters 
were  obtained  by  sampling,  additional  sampling  of  these  lead  time  demands 
may  very  well  be  warranted  to  reduce  the  associated  standard  deviations. 

Table  5 also  suggests  that  the  optimal  solution  value  is  very  sensi- 
tive to  all  of  the  item  costs.  If  the  structure  of  Problem  (SC)  is  examined, 
tHis  result  may  at  first  appear  contradictory  since  the  c.  appear  only 

| 

I 
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in  the  investment  constraint  and  the  optimal  value  function,  according  to 
Table  5,  is  apparently  not  very  sensitive  to  the  investment  limit  . 

The  problem  is  one  of  precise  interpretation.  The  partial  derivatives 
measure  rate  of  change.  But  inspection  of  the  investment  constraint 
g^(Q,r)  at  (Q*,r*)  reveals  that  the  change  in  an  item  cost  c^  by 
any  amount  Ac^  has  the  same  effect  on  the  constraint  as  a change  in  the 
investment  limit  of  -(r*  + Q£/2  - p^)Ac^  . Since  the  quantity  in 

parentheses  may  be  verified  from  Tables  3 and  4 to  be  much  greater  than 
one  for  all  i , it  follows  that  the  effect  of  changing  any  c^  by  any 
increment  6 will  be  much  greater  on  the  constraint  (and  hence,  on  the 
optimal  value  Z*  , since  the  constraint  is  binding)  than  the  effect  of 
changing  by  the  same  amount  6 . This  implies  that  | 3Z*/ ^c^ | > 

|3Z*/9K1|  for  each  i and,  in  fact,  it  can  be  shown  here  that  9Z*/9ci  = 
-(r*  + Q*/2  - i0  9Z*/9K^  , so  that  the  relationships  indicated  are  indeed 
precisely  verified. 

The  above  observations  might  also  suggest  that  some  care  must  be 
taken  in  interpreting  the  results.  Changes  in  the  parameter  associated 
with  the  largest  (in  absolute  valve)  partial  derivative  will  give  the 
greatest  local  change  in  the  optimal  value  of  the  objective  function, 
compared  to  a change  of  the  same  magnitude  in  any  other  parameter  taken 
individually.  This  follows  because  either  the  objective  function  and/or 
some  of  the  constraints  (as  above)  are  most  significantly  affected  by  this 
parameter  change  ar  the  current  solution.  General  rules  have  not  been  given 
for  selection  of  optimal  changes  in  the  parameters,  i.e.,  for  determining 
the  optimal  magnitude  and  combination  of  such  changes.  It  is  well  beyond 
the  scope  of  this  paper  to  pursue  this  "macro-analysis"  determination, 
though  it  should  be  noted  that  the  greatest  local  rate  of  decrease  in  the 
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optimal  value  function  is  along  the  direction  of  the  negative  of  the 
gradient  of  this  function  in  parameter  space  (i.e.,  along  the  vector 
composed  of  the  negative  of  the  components  of  the  partial  derivatives 
with  respect  to  the  various  parameters).  A user  would  nonetheless  have 
to  determine  the  feasibility  of  this  direction  of  change  and,  if  feasible, 
the  optimal  move  along  this  vector,  taking  into  account  other  factors  such 
as  the  relative  "cost-effectiveness"  of  any  schedule  of  changes  in  any 
model  parameter. 

Referring  back  to  Table  4,  we  note  that  the  Lagrange  multiplier 
u*  is  much  greater  than  u*  . Recalling  the  "sensitivity"  interpretation 
of  Lagrange  multipliers,  which  holds  under  the  present  conditions,  it  fol- 
lows that  u*  = -3Z*/3K^  and  u*  = -3Z*/3K2  . This  conclusion  is  consis- 
tent with  the  result  obtained  in  Table  5,  and  it  means  that  the  work  load 
constraint  g^  is  by  far  the  more  effective  in  determining  the  minimum 
number  of  expected  time-weighted  shortages  at  the  current  value  of  the 
parameters,  e.g. , a small  increase  in  will  have  a greater  effect  on 

reducing  Z*  than  a small  increase  in  . Nonetheless,  a user  must 

again  simultaneously  consider  the  comparative  costs  involved  in  making 
finite  changes,  in  conjunction  with  their  expected  effects,  to  arrive 
at  an  optimal  marginal  improvement  based  on  this  first  order  information. 

The  sensitivity  information  is  valuable,  but  requires  some  care  in  exploiting. 

Table  6 gives  the  estimates  of  the  first  derivatives  of  the  optimal 
reorder  quantities  and  reorder  points  r^  with  respect  to  each  of 

the  problem  parameters.  This  is  extremely  detailed  information  which 
gives  an  indication  of  how  the  components  of  the  solution  vector  itself 


will  change  as  the  various  parameters  change.  In  particular,  this  information 


25 


T-340 


SOLUTION  POINT  SENSITIVITY 


PARTIALS 

ITEM  i 

1 

2 

3 

3Qi/3K2 

- 47.3187 

- 18.7610 

- 14.9065 

3ri/3K2 

5.2265 

6.1961 

9.9671 

3Qi/3c1 

-208.8688 

15.3140 

14.3755 

9ri/3c1 

- 31.7918 

- 10.3425 

- 20.0020 

3Qi/3a2 

.8469 

. 2271 

.1084 

9ri/  9a2 

.1273 

1.0337 

.4919 

3V3C2 

8.1908 

- 4.9719 

3.8522 

3ri/3a2 

- 2.6783 

- 7.2676 

- 7.1087 

3Q  ^303 

- 1.2374 

.4033 

. 5843 

3ri/30j 

.2611 

.3839 

.0446 

3V3c3 

.1670 

.4442 

.4251 

3ri/  3c3 

- 2.3072 

- 3.1702 

- 12.1523 

can  be  used  to  obtain  a first  order  estimate  of  the  solution  vector  of 
a problem  involving  different  parameter  values,  having  obtained  a solution 
for  a given  set  of  parameters. 

The  partial  derivatives  of  the  Lagrange  multipliers  with  respect  to 
the  parameters  are  given  in  Table  7.  Again,  these  can  be  used  to  obtain  first 
order  estimates  of  the  Lagrange  multipliers  of  a problem  with  different  pa- 
rameter values.  In  particular,  the  relative  effects  of  the  constraints  on 
the  optimal  value  of  problems  involving  different  parameter  values  can  be 
estimated.  Furthermore,  it  can  be  shown  that  the  partial  derivatives  of  the 
multipliers  with  respect  to  and  1^  yield  the  second  partial  deriva- 
tives of  the  optimal  value  function  with  respect  to  the  parameters  and 

K2  , under  the  present  conditions.  Thus,  the  kind  of  information  given  in 
Table  7 can  be  used  to  provide  a second  order  estimate  of  the  optimal  value 
function  Z*  for  different  values  of  these  parameters. 

To  illustrate  and  test  the  application  of  the  type  of  information 
provided  here,  the  first  partial  derivatives  with  respect  to  c^  of  the 
optimal  value  function  Z*  , the  solution  components  and  r^  , 

and  the  Lagrange  multipliers  u^  and  u^  , were  used  to  give  a first  order 
(Taylor's  Series)  estimate  of  the  corresponding  solution  values  associated 
with  the  problem  where  the  given  value  of  c^  was  increased  by  one  dollar. 
These  estimates  were  compared  with  the  respective  values  of  the  solution 
obtained  by  actually  solving  the  perturbed  problem.  The  results  are 
summarized  in  Table  8.  Though  the  perturbation  is  large  (the  parameter 
being  increased  by  100%  of  its  current  value),  the  estimates  are  seen  to 
be  extremely  accurate  with  the  exception  of  the  estimated  reorder  quantity 
. Many  uses  could  be  made  of  the  estimated  solution;  e.g.,  should  it 
be  desirable  to  solve  the  perturbed  problem  accurately,  it  would  be 
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TABLE  7 

L.M.  SENSITIVITY 


OPTIMAL  L.M.  DERIVATIVES 


PARTIALS 

CONSTRAINT  i 

1:  INVESTMENT 

2:  WORKLOAD 

3ui*/8K2 

-•0002 

-.1382 

3ui*/3c1 

. 0006 

'.1635 

3ui*/3o2 

.0000 

.0006 

3ui*/3c2 

.0002 

. 0489 

3ui*/3a3 

.0000 

.0020 

3ui*/3c3 

.0003 

.0323 

computationally  extremely  advantageous  to  use  the  estimated  solution  as 
a starting  point. 

The  complete  exploitation  of  this  sensitivity  analysis  information 
now  available  will  depend  largely  on  user  interest  and  ingenuity. 
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TABLE  8 


FIRST  ORDER  ESTIMATES  FOR  A UNIT 
INCREASE  OF  PARAMETER  c. 


QUANTITY 

i 

F(0 

ESTIMATE 
u t 
F(e  )2 

ACTUAL 

' Z ABS. 
ERROR 

z*^1) 

15.159 

14.996 

1.08 

QjCe1) 

324 

412 

21.36 

r1(e1) 

221 

229 

3.49 

QjU1) 

261 

257 

1.56 

Tz^1) 

267 

268 

.37 

QjU1) 

299 

297 

.67 

r3Cel) 

417 

420 

.71 

UjCe1) 

.0058 

.0057 

1.75 

UzCe1) 

.7865 

.7671 

1.94 

TF(e1)1  = F(e°)  + (e1-e°)TVEF(e0) 


= F(e°)  + (l)aF(E0)/3C;L  . 
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APPENDIX  A 


SUBROUTINES  SENS,  LMULT  AND  PRESEN 


SUHOQoTINe.  SENS 


IS  MfiCH  19*2 


This  VERSION  or  THE  SIMM  T I VI  TV  ANALYSIS  SUMPOUtlNf  IS  US*0  TO 
COMPUTE  TK  O I c r C T l ('  N 4 L Cl  «•  I VAT  I Or  X AND  r WITH  P'lSPECT  TO 

CFfc  T A l N P*HMrTF»-5  COOt  O IN  THC  AOFAY  P*B(?f  J,  THf  OIO^CTIONAL 
0*.«I  VAT  I VCS  AG*-  EST  I MA  T c0  FJC  ON?  PAfclMF.T'iB  AT  A T I Mf  WITH  NPAP  Bf.  I NO 

The  nuw'>:c  of  pacamctf^s  involved  in  the  sensitivity  analysis*  the 

USE  OF  THr  PA«AMTT*tS  PAW(?')  MUST  9“  CONSISTENT  THC OUGHOUT  THE 
USE  p • S SUr-OUT INFS. 

Tnr  SU^COUTINE  IS  USED  FOC  A SENSITIVITY  ANALYSIS  AT  THE  FINAL  SUB- 
Pec^LtM  C*-  F 0*»  A SENSITIVITY  ANALYSIS  AT  EACH  SUPPFORLe.M  ALONG  THE 
•HM“.7|N5  TGAJECTC^Y.  po*c(2CI  TS  T»l#  AOGAY  OF  Olf  FWESClHG 
Jf4TEPV*LS  C?FU«;  SPONf^I  NO  TO  TH"  P4PAM=T‘.F$  PAP  ( 20  I • DPAPIgO)  IS 
ASSIGNS!/  VALUrS  IN  SU»*C  OUT  IN^  PAOOIF. 

This  A^r.f.ACH  TO  SENSITIVITY  ANALYSIS  IS  OUC  TO  A*  v.  FlACCO*  The 
F1PST  Vf-SION  WAS  ccpro  RY  p.  CAUSCY.  THE  SECCNO  VEFSION  WAS'COOrO 
BY  *•  C*  m VL ANOE° • THIS  IS  THE  TM|®0  VERSION  WHICH  IS  AN  EXTENSION 

GF  THC  «.*  CCNO  vrcsfEN  t r PFFMIT  SrNS|TJvITr  ANALYSIS  ALONG  THE 
MIN|M|/1>JG  TP*jt  CTCPVt  AND  «AS  COOrO  MY  P.  L.  APMACOST# 

IMPLICIT  FE*l*«( A-h.O-Z) 

FE*l*4  FhGI'4«PAt  lO»f  PSI  |THC  T*C  i KrPl  • XEP? 

CCMVC*4/SHAE-;/X  C2rJ.LrLC20».Af?3.23».N,M,  «N,NP1  * NM  1 
CCmmcn/SOAL /H  « HI  ,i»2 

CCMVCN/CPTNS/NTl *NT? ,NT 3 ,NT« .NTS . NT6 .NT7 ,NT6*NT9 «NT1 0 
COmn  >N/VALUE/F  .G.pr  , ssiG^A  , F jf  a<»  ) ,whq 

CCvmL'N/CcST/0*L*<  2F).orLXri?01*FHOlN.FATIO.rPSI.THETAO* 
INTClc.NUMJNI  • < 1 Or(*)  2:  t ,n2t?r  >#x-?CPC).xnil?0)*PPI* 

»PC?,PI  .FI  ,Oj|  (40  I ,OOTT .PGF A?|  20  ) . 0 1 AG(  20  ) . 

3Pf  £V3.»0cL*.tSIGI .Cl  • NPhASE »NSA  T I s 
CO-4***  N/SFN/O  AP  ( ?c  ) ,r>PAG  | ?;  ) .NPAP  . I Sc  NS 

CCMMjs/rXPCPT/».CF04*|  .NCX0P?.NEX0P3.NEX0P4,NEX0P5*XEPI  • XEP2 
01  PC NS  I r N DLL T I . OTLT *( 2 3 > . XFHl 20  » , X3HI2C I 
0 I +* NS  ION  DEL  M'j ( «0  > *OU(40  I «KTCST  ( 2?  ) 

»«f!lCP4  * | 

NCXOn*  > 0 

OO  S 1 * I • N 

Of  l IXC  I I ■ OLLXI  I | 

9 Of LT | | ) « DEL  *0 ( I ) 

cail  stijsf 

aval  * i.cr-io 

PPM l * * 4 4Z 
IS*. NS  » I SENS  ♦ 1 
CALL  PA4-CIX 

wg i t c out  x.  chc  and  pap. 

•BIT*  ( S • * 0 ) cmo 

|0  FCFMAMIH  //T3X.2GMSI  NSIT  I VITY  ANALYSIS  // 

IS*  • 22** T M*  VALUE  CF  PlfHQI  IS  .EI2*S) 

WPI  Y=|  A.,201 

JC  FUFNAfl/NX,  A4HT»ir  POINT  AT  which  THE  ESTIMATE  of  sensitivity  WILL 
I tiE  MAO r IS  ) 

00  4?  |t| • N • 6 

% I f *MIN*)I  I #S.M 

»®ITr|/,.j0|  (U.IUII.  Jtltll) 

33  F OF m 4 f ( 6<2X.2HX< • f ?.?H|«,G|4 .7)  f 

• C C ONH  Ml* 

VMTM6.50)  4 i I .PAF  I I I ,f.PAP(  I ) ) , HI.NB4D) 


Cop y available  \o  DDC  cioc-a  M 

jvrjs/t  hiuy  leaib^  repto ducted  , 


CO  3 J 

0014 
00  35 

CO  36 
0037 
CG36 
CC39 
0040 

C04  j 

0042 

0043 

0044 

0045 

0046 

0047 

0046 

0040 

0050 
CC  51 

0052 

0053 

0054 
CC55 
00  56 
0057 

0056 
OC  50 
CC50 
C061 
C C 62 

0063 
C C 64 
0065 
CC66 
CC  67 

CCE  9 
C 060 
CC70 
00>l 
0072 
00>3 
C07* 

0075 

0076 

C 0 77 
0076 
00»9 
0060 

0041 
0062 
0043 

0064 
0365 
00  66 
OC  67 
0364 
0059 
0093 
0091 
0072 

0093 

0094 

0095 

0096 
CC97 
0096 
0099 
0190 
0101 
0 102 
0103 
0194 
0105 


50  FO*'*aT  4 46M)  r*UMrrco  VALUE  O16FERENCIN0  INTERVAL  ✓ 

1 (|A*?K«0|4,6,6K.G14.5I  ) 

1 TL  <f  .551 
55  F C*-MATC//| 

C EVALUATE  f,  G 4N0  M. 

CALL  *•  c.  SINT  ( 3.F  | . V 

|r<MI'M/#f  0.0  1 GO  TO  A 5 ♦ 

00  43  J=1 , 

CALL  4f  STNT 4 j,Rj( J)  | 

60  ccntinuf 

IF  4 M XOI-n.CO#  31  GO  TO  05 
CALL  PPE6EN4 OU.KTEST ) 

C COHPU»r  OiL  F. 

65  CALL  GP A C 14  0 1 

on  r.  i*i,n 

X3M4 l 1 * OEL<  I 1 

99  CONTI noC 

C Cl)4POTr  COfL  1**2  P - 5TQOEO  IN  A# 

CALL  SEC0P042I 

C PERFORM  ThC  l-U  DECOMPOSITION  OF  A # 
on  ir 3 t«l  ,N 
Oc  L * 4 I I *C  • 9 

loo  ccntinuc 

NP * N T 3 
N T 3 * 1 

CALL  1NVCCS4 1 I 

C CMFC*  TO  *AAKE  SUP*  AN  ORTHOGONAL  MQvC  IS  NOT  ATTEMPTED# 

00  I 1 C I * I .N 

IFCD*  LX4  1 I.EO.O.CI  GC  TO  HO 
VRITC 46 .1051 

1C5  rCPMAT494Mi  THf  matrix  OF  SECOND  PAPTIALS  IS  NOT  POSITIVE  OEFINl.TE 
I.  Sc*»  S I T t V I T V ANALYSIS  IS  TERMINATED  I 

GO  TO  ?'? 

11C  CONTINUE 

DO  l«:n  1*1  »NPAC 
X2H4  I I » PAC  4 11 

123  CCNTJNU'  f 

OO  ?f  9 Jr  1 ,NPAP 
|F  4M  *LP*  .’I*  .4  | GC  TO  116 
IF  f K*r  si  ( J I U.3  ) Of)  TO  2 00 
115  IF4K  X.'PA.t'J.JI  GC  TO  121 
CALL  L^ULTI  7.rCL*fU.0*:M,f.O» 

C COMPUTE  0(C»CL  PI/DMJ)  AND  0(F)/04(j)  USING  CENTRAL  O I FFE  PENC 1 NG< 

121  PAPTjl  *PA=(J»  ♦ OP  AP  4 J 1 

CALL  P r S T N T ( 3, OF  ) 

IF(M#E0#C)  GC  TO  126 
00  126  I s j <*  **  • 

CALL  Pr  STNT (l.Pj(l)) 

IF|C J4  I l.GT. AVAL  I GO  TO  12* 

123  OPAP  ( J)  «f  . 1 *r.PAR  4 J» 

PAR ( J ) * X2H4JI 

«PITf(e«124|  J.OPAOfJ) 

124  FOPWAT416H  PCSFTTING  OP»vP.  (i  I\2v3h>  ■ • G1  4 • 5 ) 

IFCOPAkCJ) • E 0 • 0 • 1 GO  TO  291 

• GO  TO  121 

125  CONTINUE 

126  IF4M2.EQ.O)  GO  TO  126 
OO  127 

M/P| «M4 | 

CALL  RESTNT|*/p| .PJ4M/PI I I 

127  CONTINUE 

126  IF  4 N«  «0M4#E0#c.  I GO  TO  129 
CALL  L*ULT 4 1 .0EL*U.0EM.0U> 

129  CA1L  GMC42I 
OO  l?0  l*t*N 
DLL X< I I aOFLXOC I 1 
ISO  CONT|NUr 

Of  M*2#t APPAW  < j ) 

P*M  Jl  ( J 1 - OEM 

CALL  PTSTNTtO.VAL > 

IFIM.r.O.ri  GO  to  136 
OO  175  I * I • M 

CALL  t r S T N f ( I • P J 4 III  # 

If |OjC  I I . GT. AVAL  I GO  TO  135 
GO  TO  12 J 

135  CONTINUE 

136  |F|M/#?O.0l  GO  TO  134 
OO  137  1*1. M2 

M/P| *M4 | 

CALL  RE  SINT  4 m/p|  .PjfM/PI I I 

137  CPNTJNUr 


0 172  30 
01 7240 
C 17250 
01 72*0 
017270 
01 729C 
01 7290 
01 73C0 
017310 
Cl  7 320 


017330 
01 7340 
01 7350 
Cl  7360 
01 7370 
01 7340 
017399 
01 7490 
017410 
017420 
017430 
Cl  7440 
017450 
Cl  74  60 
Cl  74-»0 
01  74*0 
017490 
Cl  7530 
C 17510 
C 17579 
017530 
017540 
017550 
01 7560 
C I 7570 
017550 


Cl  7590 
Cl 76C0 
017610 
Cl 7C20 
Cl  703C 
Cl  7f  4 0 

Cl  7650 
Ol 7660 
01 7670 
017040 
01 7690 
C 1 7700 
017710 
017720 
©17730 
017740 
017750 
017760 
C 1 77  70 


017790 
01*600 
C I 74 1 0 
017620 
Cl  • e «v 
01 7640 
Cl  7650 
Cl  7 66 C 
01  76  T 0 
01 *659 
017090 
01 7930 
017910 
01 7920 
01 793C 
017940 
01 7950 
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nkitd  II# 


134 


0106 

C107 

0 l 34 

01  10 
0111 
0112 
0110 


0114 

01  15 
01  16 
0117 
0114 
0116 
0120 
< 1 21 
0122 
0121 
Cl  24 
Cl  25 
C I 26 
0127 
0120 
f : 29 
01  10 
01  31 


l 39 


MO 


IFI  N'.HCf  AMQ.CI  CC  TO  139 
COLL  L"OLTC  2 .DtLMU.OFM.PU) 

CALL  G*-*0<?1 
C»F  * f DF  - V AL  )/0*M 
00  MO  1*1. N 

Oft  XI  1 M105LX1  1 I - DELX01  1 1 I/OEM 
CCNTlVOF  » 


OAF  ( J I «•  Jt?H|  J I 

C HAVING  FACIOMD  4.  SOLVE  A*X*8  FOR  X* 

C KHiftf  H-O‘1*  »N0  X?DX/DA|J>. 

CALL  INVMSl*) 

C PF  INT  OUT  t 4/0  A t J I 

• F f T/  </.  , | 5f  | j 

150  FQ6MATI47M  X-Or F 1 VAT  I Vf 5 ABE  WITH  RESPECT  TO  PARAMETER  .12) 

00  1 70  1=1 ,N,6 

■ MU  I'  ,|'0)  C J J.DfLM  J J>  . ^ 

160  FOFV*T4  « (ah  0X4 • 12 */M| **G1 4.7 ) 1 % 

I7C  CONTINUE 

IFlNf  XCTA  .*-0.0  ) C.O  TO  375 

ir lY.rt.r  ) oc  ic  a."'! 

CALL  LA'ULfl  ! .O^L^U.Cf  m,pu> 

■*11*16.350)  J 

350  FORMAT ( ✓ AIM  U-OFR | VAT | VF S KITH  RESPECT  TO  PARAMETER  •(?! 

OO  T* 1 I * I ,m  ,6  / 

11  * MJNCC  M5.*»l 

• R I Te  I / .752)  U JJ.Oui 

352  rOFMAMtMH  OUl  • 1 2 . 2H  ) * • 0 ) A • 7 ) > 

351  CONTlNUF 


017070 

017940 

017990 

oncoo 

oiecio 

Cl  AC  20 
01 6030 
ci s: ao 
016050 
OM  60 
01  PC  /o 
oi eoeo 

01 90  90 

0191  CO 
019110 
0161 20 


0 1 32 

331 

|F4*£.FQ.C)  CO  TC  375  • 

01  33 

CALL  L MUL  T 1 4*  .0*  L *U  • Ct  M • CO  ) 

61  34 

■ F I T«-  4 6, 360  ) J 

0125 

360 

FUFMAT4/  4lr  ■ - OF  P l VA  T 1 vE  S WITH  RESPECT  TO  PARAMETER 

C 1 '6 

DO  361  Ml. M2, 6 

0137 

I 1 • MINCI  MS. M2) 

! 


Cl  34 

■ RITMN.267)  1 1 JJ.OUl JU*m) > • jjr[,n  ) 

CI  39 

362 

fOFMAT46(AH  owl  • 12..2HM.G14.7)  ) 

CMC 

361 

CUNT INU5 

0141 

375 

COM  1 nu« 

C COMPUTE  OF /DA f J ) • 

C19I 3D 

C 142 

00  in  Ml.N 

C I PMC 

0143 

or  * OF  ♦ x )Hl I ) *DF L xi 2 I 

C l 4 1 5 3 

0144 

1 6 C 

CONTINUE 

C14163 

c pf  i»it  or/CA  ( 5) . 

Cl  41  7 D 

CMS 

«* 

«R|  TP  4 6. 1.90  I OF 

C 1 31 4? 

3 1 A 6 

19? 

01  41 RC 

Cl  4* 

200 

CCM  1MLE 

C 1 92  DC 

Cl  44 

GO  TO  Pr  2 

Cl  4210 

C 1 4 9 

201 

• P1TF1F , 1 Co  ) J 

01 822C 

CISC 

106 

FORMA  T 1 |H  . 21  MTF.RM1  NAT  I NG  P AR  4M*  T ER  . 1 3 . 1 6H  OUE  TO  0®AP  ■ 0 

/) 

019230 

3151 

GO  TO  200 

C 1 62*3 

C 152 

202 

N1  J*NP 

C162  5C 

0153 

CALL  REJECT 

01*260 

C 1 54 

00  205  1* 1 .X 

C 1 P270 

Cl  55 

OGLXl  1 I a 0ELTX4  1 ) 

C 1 6293 

0156 

205 

OtLX Y<  l ) * OEL  T < 1 1 

Cl  4253 

0157 

PF  TURN 

01 9330 

0159 

€ NO 

016310 

| 
I i 

- 


1 


/ 


1 

i 

I 

i~ 
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SUt'OClUl  |wf  L*ULT|  ! KD  . O*  L mu  • OC  m . OU  I 

ILLICIT  P*:  *L  •*»(  A-H.O-V  ) 

Cfif*  ft  hC ! N, -*  T 10.  rrs  I .THE  T4C  .KfPl  . Jtr^2 

COMMON/SHAC ~/K( PC  I. OF I 

CC^M'IN/VAL  UC /r  , G.Pfl  • C S I r.MA  , C J ( 4 0 I ,«>MQ 
CO^M.iN/fTOAL/  H,  HI,  »l 

CCMMUM/CP ST/OrL* ( 20 ) .Of  LX' ( ?Q)»PWCIN,OATIO*FPSI  ,T  HE  T AO 
iNTCTft  , KU^IM  • *1  < 20  • X2(  2?  I . X 3(?0  > »*F  2«  *0  » • XF  I (20  I «PH|« 
2PP2  ,p  1 ,F  1 ,0  J1  c 4 0 I .OPT  T , Pc.C  *0  ( 2?  ) .0UC(?0), 

3PPf  V?  . AOELX  .«SI  Cl  • r,  1 , NPH*SE  , NS*  T I S 
OI^^NSION  OKLMUC«0)  .liU<A*l 

MM2  s M * M£  ^ k 

GO  TO  (1,2.3.11.  1ND  ^ 

c t ho  * r 

OO  too  I * I • MM  2 
iro  ort m«j«  1 1 » o. 

PETUftN 

C 1*0  « I * 

1 OO  50  I « 1 .«*«/ 

50  ortMiii  ii  c ft  j(  I » 

PC  TUAN 
C INO  * 2 

2 DO  ^ !•!, “Ml 

60  Ofl^U(l)  * lOTLMOCII  - BJ(|))/OEM 

PLTU=N 

3 OO  Tr  | «|  , mm 2 
CALL  G ft  A 0 I ( I I 
CALL  PLSTMT(  I , VAL  ) 

SUM  e 0. 

OO  71  JJ*|.N 

71  SUM  * SUM  ♦ OcL< JJIAOELXC JJ) 

lf< IND.EO.A ) GC  TC  PZ 

OU<l)  *— ( SUM  ♦ 0-LMU(  I I ) *ftHO/( V*L*«2 ) 

GO  TO  70 

60  OUll)  ■ ?,*CSUM  ♦ DFLMUC I I »/ftHO 

7 0 COOl  HO£ 

BE  Tuft H 
E MO 


COOl 

COO? 

COO) 

OCj* 

0005 

0006 
0007 


0 0 0* 
C030 
OCIO 


001  I 
0012 
0013 


0016 

0015 

0016 


C01  7 

oc  la 
cc  I 9 
OC  20 
OC  21 
CO  2 2 
0 023 
0C?a 
0?  25 
CC  ?6 
CC  ? 7 
cc  ?e 

00  29 
00  30 
00)1 


CC  23 

0 C 24 

30 

0G25 
CO  26 
OC  27 

40 

CG20 

50 

0029 
OC  30 
CO  31 

60 

CC  32 

?C 

CO  33 
0 C 34 
00  35 
C 0 36 

60 

00  37 
CC  36 

90 

0 0 39 
C040 

004| 

0042 

95 

CC  4 3 
0044 

100 

0C45 
OC  46 
0047 
0046 

600 

0C49 

601 

005C 

0051 

005? 

CC5J 

00  54 

0055 

200 

CC  56 

005? 

0056 

2 50 

C 059 
C 0 60 

602 

OC  51 
CO  62 

603 

0063 

0064 

604 

0065 

300 

0066 

0C6? 

605 

nnrV».H'U,Kl»sr) 

IMPLICIT  Ci  4L*S<4-m,C-ZI 

2f  I *r»rLC  23  I • M*V  i ?C)  .N.m.mn.NPI  ,nmi 

cr**  ,N/V*L  Ur/F  PC  .*  S I C.MA  J4  40  ) .CMO 
COMM:i-</f  Qk[_  / M,  Hi  «M£ 

CC**CN/cr  N/o*o  4 20  > . 0M4«-  120  1 * NPAR  , (SENS 
01  MCNSIC,r«  G«<ac)  .OU(40)  .KTEST(20).KLlST(20) 

MPW/  * W » M / 

FT«$!  * C .f  J!  • CAHSIF) 

00  l?0  J = I . sp»; 

<TF  ST ( J ) * J 

PAR|J)  s o*t,  (j)  • OPAR(J) 

call  er.STur  ( o .or  i 

ir  IMP-4/.  r&.  , CO  to  20 

00  1 * 1*1  .•«.»*£ 

CALL  «rSTNT | I ,DU<  I I > 

Continue  , 

OEM  * 2.  • 0P*MJ> 

P*M  Jl  * o«o<j|  - OCA 
CALL  «-r5TNT(0.*r) 

ir <»p»z,£g.i j cc  to  *o 

DO  33  1*1  .mpm£ 

CALL  cr  5 T N M f .0x4  lit 

costinu* 

or aps  - cor  - xfi/pfm  ' 

ir(Aav2,ro.3)  GO  TO  60 

00  *50  1*1  ,«P"l 

oucii  = coucri  - oMm/ocA 
su*  * offps 
if4w.eo.oi  gc  to  eo 
oc  ?3  i*i  .* 

sum  * suM  - rmc/pjc i ) *du( i i 
I F 4 **£  • CO  • 0 ) GO  to  95 
TSOM  = 0. 

OH  9^  1*1  ,«z 

|M  r | 4M 

T SU*  * T SUM  ♦ F J( |M|*OU<  !M| 

SU*  ■&  Su*  4 ISO*  • 2./RNO 
OFLIJ)  * SUM  - 

PA*3(J)  * PAP(J)  ♦ OP  AO  C J | 

OT*ST  * C>  A 3S  ( Of  Lljl) 

ir 40TF ST .Gr .F TEST > KTEST4J)  * 1 

CONTINUE 

• R I T»  It  ) 

FORMAT  ( ?z%  . 1AH-DPT  IMAL  VALUE  FUNCTION  SENSITIVITY  //I 

DO  20  C 1*1  .NPAh  ,5 

1 !*MIN?|  | *4  . KOAC  | 

• 0IT-C6.FC1)  < ( JJ.0EL4 JJ»  I . JJ*  I. II) 

roKMATC  547M  OT/OA( . I2,2m)e,gi*,7|) 

CUNT  INU1 

JJ  * 0 ' 

00  250  J*1,NPAP  ' S 

IFCKTC STC J) .FO.C I GO  TO  250 
jj  * jj  * t 

XL  I ST  C j J ) • J 

continue 

IFCJJ.r.0.0)  CO  TO  300 
•PITCC6.602) 

FORMATI  /5IH  CETAItfO  SENSITIVITY  RESULTS  FOLLOW  FOR  PARAMETERS 
WRITE (6.6C3)  (KLISTCII*  1*1. JJ) 

FULMATIJM  40(12. ?h  •)) 

• »l  T*.  C t .<  04  ) 

FORMAT  < / ) 

FT  TURN 
WRITE 46.105) 

FORMAT!  4.*m  ThCPF  A*»R  NO  DFTAIlCO  SENSITIVITY  RESULTS  /✓! 

RETURN 
NO 
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APPENDIX  B 


USER  SUBROUTINES  FOR  SCHRADY-CHOE  PROBLEM 


©OCI 

SUBROUTINE  CCA01N 

CC  32 

fMfLicrr  cf4L*«(  a-n.o-z  » 

0003 

Cn*'*'ri*i/  l sy/'lG  T ^ < pq  i ,i>h|  < ?n  | ,i>rN5E(20  ) • I OCN*  ( 20  1 * Nl 

0C3* 

COMMON/ SEN/P A3 < 2C I • DP AP ( 20  1 .NPAR * I SENS 

0035 

901 

FCPMAT ( I5.4F12.C 1 

0006 

0E4DCS.9CII  M *PAK< 1 I .PAP <2» 

0 00  T 

■ C ITT (A.95I  ) N I , PAP ( 1 ) , PAQ ( 2 » 

0006 

00  MO  1 *1  .Nl 

0009 

FEZ  DC  5. 50  1 ) IOfNT<*I  )•  (PAC  ( A*!-2*  . J=1  .A  ) 

00  10 

*R  1 Tc  ( 6 • 90  I 1 t Of  NT  < 1 1 • (PAP  ( 4 * 1 -2*./ ) * J*1  .*> 

0011 

100 

CONTINUE 

0012 

NRAF  9 4 • N I ♦ 2 

0013 

PETUFN 

COM 

END 

ccot 

SUk'PCJUT  I NE  PFSTMI  IN*  V*L> 

0002 
C 003 

IMPLICIT  Pf AL«P< A-M.O-Z) 

C03* 

CO*«ON/SM«nc/x(;>r  ) t OF  L ( 2-3  > . *(  20*  PC  ) • N,  m,  “N,NPJ  » NM 1 

C0G5 

COMMON/ IN V/uETAC  ?0»  .PHM201  *DLNSE<20  1.  IOENT(20>*NI 

0006 

CCMVON/SEN/PAR  ( 20  > . OP  AP  ( 20  ) • NPAP  * ISENS 

COOT 

VAL  * C. 

coca 

|F( |N.f 0.0)  CO  TO  300 

OCOP 

IF ( IN, 50.1 ) GC  TO  TOO 

C010 

2C0 

00  250  1*1. Nl 

con 

U * 2*1-1 

0012 

QO  * K< I J) 

0013 

I F f 00 , L E • 9 • 1 GO  TO  180 

CC14 

250 

VAL  * VAL  ♦ PAP(**1*2»/Q0 

ec  is 

VAL  * PAP  { ? t - VAL 

CO  1 6 

RETURN 

0017 

1 BO 

VAL  *-1.0 

C 0 1 6 

F£  TUP  N 

0019 

100 

OO  ISO  1*1 *NI 

0020 

IJ1  * 2*1 

0021 

IJ  * 131  - 1 

0022 

RF  * X ( 1 J 1 ) . 

©023 

1 F(PP ,LE .0. 1 GO  TO  160 

002* 

00  * K( IJ ) 

C C 25 

if(oq«lt*o*>  go  tc  leo 

0 026 

ISO 

VAL  * VAL  ♦ PAR|**1*| )•( P P ♦ 00/2 •“ PAR «*• 1 - 1 >) 

0027 

VAL  * PARC  1 1 - VAL 

0026 

RF  TURN 

C 0 29 

390 

OO  3r 0 1*1. Nl 

©0  30 

1J1  * 2*1 

0031 

IJ  ■ Ul  - 1 

0032 

00  ■ X( I J) 

C 0 39 
00  34 

PC  « Ki I Jl ) 

0035 

UU  • RAP (*• I-l 1 

00  16 

OC  3 7 
fit* 

SS  * PAF (**I 9 

OC  99 

oelta  * CP  - UU 

00*0 

IN  • delta  / ss 

«l<  41 

Call  anotrc /n,f i ,ueni 

0 C * 2 

PM I III  * FI 

60*1 

DLNSF III  * OF N 

C9«« 

0FTAII)  * 0 • 5 • ( ( SS*SS*OELTA*OELT A 1 *PM | f 1 l-SS*OELT  A*OENSE ( 1 » » 

09*5 

♦ 350 

VAL  * VAL  ♦ OTTAI 11/00 

C0*« 

RE  TUPS 

00*7 

END 

I 
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0031 

SUH**PVJ  1 1 \r  GC*OM|N) 

c::2 

IMPLICIT  fcCAL»H< A-M.n-Z 1 

cto  j 

CC*M  3N/5M * = r:/*<  2f  ) ,(>FL<  ?9  1 • 41  20  • 20  > .N.M.  MN.NPt  .Nttl 

0C94 

COMMON/  | NV/  »•_  1 A(  70  1 .PH  | < ?0  } .DENSE  C 7 0 1 • I0FNT1  20  1 , Nl 

0005 

ec»*w-)N/«,f  N/P4W  ( ?r>  1.0P4P  120  1 •NPAfc  • ISt  NS 

C 0 06 

in iN.ro. e i oo  tc  300 

0307 

|E< 1N.F0. 1 » CO  TO  100 

C396 

00  250  1 S 1 .Nl 

CC39 

UI  * 2*1 

0C1  3 

f J « 1 Jl  - 1 1 

001  l 

oo  * * < 1 j ) • , 

0012 

OF  L 1 l J 1 1 * 0. 

er  io 

250 

D'UIJ)  * pau ( 4 • I ♦ ? > /OQ/OQ 

0014 

PPTU4*4 

0015 

100 

OO  153  1=1 .Nl 

0016 

IJI  * 2*1 

0017 

U > IJI-  1 

CC1  * 

OCt  1 I Jl  I » -P6B  c 4 • I 4 1 ) 

oc  to 

150 

DEL  1 IJI  s OfcLl  IJI  )/2.  • 

0020 

P*  TUCN 

0321 

30  0 

OO  353  1*1. NI 

CO  22 

IJI  * 2*1  % 

OC23 

IJ  * IJI  - 1 

0024 

00  * Hi i J» 

C 0 25 
C 0 26 

fp  * HI IJI I 

0027 
C 026 

UU  • PAP(4*1— 1 | 

» 

0C29 

SS  * P AP ( 4 • l 1 

OC  33 

0331 

OFLTA  * BP  - uu 

CO  32 

IN  * DELTA  / SS 

C 0 33 

CALL  ANPTCf /N.F | ,Df N» 

00  34 

PHI (11  « PI 

C 0 35 

0!  NSF  <11  * DEN 

0336 

OE  T A | 1 ) * 0 • 5 * ( < SS«SSA0ELTA40£LT a ) 4PH! < 1 )-SS40rLTA*0CNSE< 1 1 > 

OC  37 

osl (iJi  > = <orLT4»PHi 1 1 >-ss»orNsr< i > >/oo 

CO  38 

350 

OCL(IJ)  * -or T A( 1 ) /QO/00 

0030 

•fTUAN 

0040 

ENO 

0CC1 

S060p^T|A.r  *4ATP>  I 1 ( 1 N.  |KK| 

CO  3 2 

IMPLICIT  P‘4l*A|a-h,0-2) 

CC  33 

CQvm  >n/SmA4F/<  ( ) ,CTL<  20  1 . AC  23  • 23  1 . N.M.  MN.NPl  .NM1 

00  34 

COMMON/  |NV/JeT*(?oi,PMH?CI  * DENSE  ( ?0  1 • I DENT  ( 20  ) . N! 

0305 

COAAOV/S  rn/DAf.  (2C  ) (2L  ) .NPAP,  I SENS 

00  36 

If  < 14.L0.9I  CO  TO  300 

0037 

ir  < lN.ro. i l CO  TC  i oo 

00  33 

OO  253  1*1. NI 

re  30 

IJ  • 2*1-1 

0013 

OO  * *1 I Jl  • 

oon 

250 

AllJ.tJt  * - 2 • 4PAP < 4 • I A2 1/00443 

OC  1 2 

PAT  ur-  N 0 

0013 

130 

IKK  * 1 

00  14 

CL  TUPN 

C01  5 

300 

00  333  1=1. Nl 

OC  16 

IJI  = 241 

0017 

IJ  * IJI  - 1 

0016 

OO  * X < 1 Jl 

OC  to 

PC  ■ *< 1 Jl 1 

0020 

* 

C 0 2 1 
C 3 22 

uu  4 PIC (4*1*1 ) 

C C 2 3 
OC  24 

S3  ■ PAP (4*  I 1 

CO  25 

DELTA  * CC  - uu 

OC  >6 

IN  * 0‘LTA  / ss 

C 3 2 7 

CALL  AM)TP<  2N.F  | ,rr.N| 

007*1 

PHI < 1 | r FI 

C C 29 

OENSf ( 1 > * DEK 

OC  33 

OC  T A<  1 1 *0.54<<SS4SS*OFLT'A4  0*LTA14PH|<|  | - SS*DFl  T A 4 DENSE  ( 1 1 1 

CO  31 

AllJ.tJt  * 2. 4|*f  T A<  1 I/GO**  3 

03  32 

AllJ.IJl  I*-(OLLTa  4 PM lit  | - SS  *0£N>E < 1 1 1/00/00 

0933 

350 

At  1 Jl . 1 Jl ) - PHI < ( l/oo 

0034 

or  TUCN 

0035 

END 

C 9 3 1 

SWNPOOT  I Nf  ANOfc  ( ITK.PMl  .DENSE  I 

5092 

IMPLICIT  PI  AL4PC a-m.C-2 ) 

OC  3 3 

At  4 O a 1 5 < * t ) 

0034 

* 

T * 1.3/11  .J*  *9 .2316  4*1  944*| 

00C3 

OENSf  * C • Jl«?  42  2 *01  «P  ( - v * • * K/2  . ) 

0036 

PM|  * DLNSr  *T*(  ( C I1..33C2744T  - l . 821  256  1 4 T ♦ 1 . 761  4 78  1 «T 
| - C • 35<  56  30 ) • T ♦ 0.31938151 

0097 

|F<  *X»  1.2.2 

0308 

I 

PH | « 1.0  - PHI 

0939 

2 

PfTUCN  0 

COIO 

CNO 

- 38  - 
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